Spatial Scripts

packrat::init("~/Desktop/Heloderma Spatial/Heloderma Spatial")

library(adehabitatHR) #for home range calculations
library(data.table) #manipulate S3 and S4 data tables

data.table 1.11.8  Latest news: r-datatable.com

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, first, last

The following object is masked from ‘package:purrr’:

    transpose

The following object is masked from ‘package:raster’:

    shift

library(ggplot2) #for graphic output
library(ggfortify) #to allow ggplot2 to read spatial data
library(grid) #to add annotations to the output
# library(OpenStreetMap) #for obtaining raster images
library(pbapply) #needed for progress bar

package ‘pbapply’ was built under R version 3.5.2

library(plotly) #for interactive xy plot

Attaching package: ‘plotly’

The following objects are masked from ‘package:plyr’:

    arrange, mutate, rename, summarise

The following object is masked from ‘package:ggmap’:

    wind

The following object is masked from ‘package:ordinal’:

    slice

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:MASS’:

    select

The following object is masked from ‘package:raster’:

    select

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(rgdal) #for converting spatial data

package ‘rgdal’ was built under R version 3.5.2rgdal: version: 1.4-3, (SVN revision 828)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
 Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
 Linking to sp version: 1.3-1 

library(sp) #for converting spatial data
library(rgeos)

rgeos version: 0.4-2, (SVN revision 581)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 
 Linking to sp version: 1.3-1 
 Polygon checking: TRUE 

# library(raster)
library(mapview)

Overall individual yearly home ranges for non/subsidized populations

Gila monster locations for all tracked lizards across Stone Canyon

myMap

768x1280 terrain map image from Stamen Maps. 
See ?ggmap to plot it.

OVERALL YEARLY ANALYSES

Plot of 100% MCP HRs against number of relocations

year <- read_csv("GM_Consolidated_ByYear.csv")

Missing column names filled in: 'X13' [13]Parsed with column specification:
cols(
  Year = col_double(),
  Gila = col_character(),
  Sex = col_character(),
  Environment = col_character(),
  Home_Range_100mcp = col_double(),
  N100 = col_double(),
  Home_Range_95mcp = col_double(),
  N95 = col_double(),
  Home_Range_95kde = col_double(),
  N = col_double(),
  Home_Range_50kde = col_double(),
  N50 = col_double(),
  X13 = col_logical()
)

# quick plot
# Graph1<-ggplot(year,aes(x=N100,y=Home_Range_100mcp,group=Environment))+
Graph1<-ggplot(year,aes(x=N100,y=Home_Range_100mcp))+
  geom_point(aes(shape = factor(Environment)), size = 4)+
  scale_shape_manual(values=c(16, 2))+
  geom_smooth(aes(linetype=Environment),colour="black", method="lm") +
  # scale_colour_manual(values=c(subsidized="cyan3",nonsubsidized="indian red1"))+
  # labs(title = "100% MCP Home Ranges")+
  xlab("Number of Relocations")+
  ylab("100% MCP Area (ha)")+
  # labs(caption = "Figure 3 | Non-Subsidized (Owl Head Buttes) vs. Subsidized (Stone Canyon) population 100% MCPs against number \n of fixes of the complete data set.")+
  theme(plot.caption = element_text(hjust = 0,lineheight = 0.9))
  # theme_bw()
Graph1<-Graph1+theme(axis.title=element_text(size = 18))
# legend at top-left, inside the plot
SCOH.hr.fig<-Graph1 + theme(legend.title = element_blank(),
               legend.text = element_text(size = 14),
               legend.justification=c(0,1),
               legend.position=c(0.05, 0.95),
               legend.background = element_blank(),
               legend.key = element_blank(),
               legend.box.background = element_rect(colour = "black")) +
               scale_shape_discrete(name  ="",
                          breaks=c("nonsubsidized", "subsidized"),
                          labels=c("Nonsubsidized", "Subsidized")) +
                            scale_linetype_discrete(name  ="",
                          breaks=c("nonsubsidized", "subsidized"),
                          labels=c("Nonsubsidized", "Subsidized"))

Scale for 'shape' is already present. Adding another scale for 'shape', which
will replace the existing scale.

SCOH.hr.fig

# dir.create("outputs") # create a new folder to hold the output files
# ggsave("outputs/SC_OHB_plot.pdf")

Plot of 95% KDEs against relocations

Overall combined 100% MCP means averaged across sex

library(Rmisc)
Means <- summarySE(year, measurevar="Home_Range_100mcp",
                          groupvars=c("Environment"),na.rm = TRUE)
kable(Means, format = "pandoc", caption = 'Overall combined 100% MCP means averaged across sex')

Overall combined 100% MCP means averaged across sex

Environment	N	Home_Range_100mcp	sd	se	ci
nonsubsidized	26	33.44231	20.518658	4.0240400	8.287665
subsidized	53	10.40151	6.948743	0.9544832	1.915311

Overall combined 95% MCP means averaged across sex

Means.95mcp <- summarySE(year, measurevar="Home_Range_95mcp",
                          groupvars=c("Environment"),na.rm = TRUE)
Means.95mcp

Set projection for mapping

CRS.SC<-CRS("+proj=utm +zone=12 +ellps=WGS84 +units=m +no_defs")

Function for MCP analysis

Function of MCP polygons used for mapping

Function of KDE analysis

kde_analysis.href.plot <- function(filename, percentage){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$EASTING)
  y <- as.data.frame(data$NORTHING)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
  ver <- getverticeshr(kde, percentage)
  area <- as.data.frame(round(ver$area,4))
  .rowNamesDF(area, make.names=TRUE) <- data$LIZARDNUMBER
  write.table(area,file="KDE_Hectares.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  kde.points <- cbind((data.frame(data.proj@coords)),data$LIZARDNUMBER)
  colnames(kde.points) <- c("x","y","lizardnumber")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.9,  y=0.95,
                     gp=gpar(fontface=4, cex=0.9), draw = FALSE)
  kde.plot <- ggplot() +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.5) +
    geom_point(data=kde.points, aes(x=x, y=y)) + theme_bw() +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$lizardnumber) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
    annotation_custom(units)
  kde.plot
}

Function of KDE polygons for mapping

kde_analysis.href.polygon <- function(filename, percentage){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$EASTING)
  y <- as.data.frame(data$NORTHING)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
  ver <- getverticeshr(kde, percentage)
  ver@proj4string<-CRS.SC
  area <- as.data.frame(round(ver$area,4))
  .rowNamesDF(area, make.names=TRUE) <- data$YEAR
  write.table(area,file="KDE_Hectares.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  kde.points <- cbind((data.frame(data.proj@coords)),data$YEAR)
  colnames(kde.points) <- c("x","y","year")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.9,  y=0.95,
                     gp=gpar(fontface=4, cex=0.9), draw = FALSE)
  ver
}

Function of raster of UD

# kde_analysis.href.raster <- function(filename){
#   data <- read.csv(file = filename)
#   x <- as.data.frame(data$EASTING)
#   y <- as.data.frame(data$NORTHING)
#   xy <- c(x,y)
#   data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS.SC)
#   xy <- SpatialPoints(data.proj@coords)
#   kde<-kernelUD(xy, h="href", kern="bivnorm", grid=1000)
#   kde<-as(kde, "SpatialGridDataFrame")
#   kde@proj4string<- CRS.SC
#   kde
# }

Function of trajectory analysis and distance over time

traj_analysis <- function(filename){
  relocs_data <- read.csv(file = filename)
  relocs <- as.ltraj(cbind(relocs_data$EASTING, relocs_data$NORTHING),id=relocs_data$LIZARDNUMBER, typeII = FALSE, date=NULL)
  relocs.df <- ld(relocs)
  relocs_dist <- as.data.frame(sum(sapply(relocs.df$dist, sum, na.rm=TRUE)))
  colnames(relocs_dist) <- "Total Distance"
  name <- relocs.df$id[1]
  row.names(relocs_dist) <- name
  relocs_units <- grid.text(paste(round(relocs_dist,2),"m"), x=0.9, y=0.9,
                            gp=gpar(fontface=3, col="black", cex=0.9), draw = FALSE)
  reloc.plot <- ggplot() + theme_classic() + geom_path(data=relocs.df, aes(x=x,y=y), linetype = "dashed", colour = "red",
                                                       arrow = arrow(length=unit(.5,"cm"), angle = 20, ends="last", type = "closed")) +
    geom_point(data=relocs.df, aes(x=x, y=y)) + geom_point(data=relocs.df, aes(x=x[1],
                                                                               y=y[1]), size = 3, color = "darkgreen", pch=0) +
    labs(x="Easting (m)", y="Northing (m)", title=relocs.df$id[1]) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
    annotation_custom(relocs_units)
  reloc.plot
}

Function of distance of time

dist_analysis <- function(filename){
  relocs_data <- read.csv(file = filename)
  relocs <- as.ltraj(cbind(relocs_data$EASTING, relocs_data$NORTHING),id=relocs_data$LIZARDNUMBER, typeII = FALSE, date=NULL)
  relocs.df <- ld(relocs)
  relocs_dist <- as.data.frame(sum(sapply(relocs.df$dist, sum, na.rm=TRUE)))
  colnames(relocs_dist) <- "Total Distance"
  name <- relocs.df$id[1]
  row.names(relocs_dist) <- name
  write.table(relocs_dist,file="reloc_dist.csv",
              append=TRUE,sep=",", col.names=FALSE, row.names=TRUE)
  dist.plot
}

Map of yearly HR shifts of a subset of Gila Monsters. Includes running MCP polygons, Fortify mcp polygons for ggplot2 by YEAR

Raw group 100% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex

Table 1 | Raw group 100% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex.

Environment	Sex	N	Home_Range_100mcp	sd	se	ci
nonsubsidized	female	11	22.063636	12.287414	3.704795	8.254797
nonsubsidized	male	14	43.235714	21.672372	5.792185	12.513255
subsidized	female	38	9.839474	6.889003	1.117544	2.264359
subsidized	male	15	11.825333	7.133668	1.841905	3.950494

Raw group 95% MCP home range means of Stone Canyon and Owl Head Buttes. Grouped by environment and sex

YR_GRP_Means95 <- summarySE(year, measurevar="Home_Range_95mcp",
                            groupvars=c("Environment","Sex"),na.rm = TRUE)
kable(YR_GRP_Means95, format = "pandoc", caption = 'Table 2 | Raw group 95% MCP home range means of raw data of Stone Canyon and Owl Head Buttes. Grouped by environment and sex.')

Table 2 | Raw group 95% MCP home range means of raw data of Stone Canyon and Owl Head Buttes. Grouped by environment and sex.

Environment	Sex	N	Home_Range_95mcp	sd	se	ci
nonsubsidized	female	6	20.600000	6.286493	2.566450	6.597270
nonsubsidized	male	9	38.988889	15.815139	5.271713	12.156592
subsidized	female	38	7.132895	4.280606	0.694406	1.407000
subsidized	male	15	9.195333	5.246499	1.354640	2.905415

RM-ANOVA for 100% MCP analyses between the subsidized and non-subsidized

summary(RMmod.year)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: Home_Range_100mcp ~ Environment + Year + Sex + N100 + Environment *  
    Sex + (1 | Gila)
   Data: year

REML criterion at convergence: 578.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51737 -0.36103 -0.04818  0.24140  3.14258 

Random effects:
 Groups   Name        Variance Std.Dev.
 Gila     (Intercept) 37.23    6.101   
 Residual             85.97    9.272   
Number of obs: 79, groups:  Gila, 31

Fixed effects:
                                Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                   -991.17075 1749.77011   67.87392  -0.566  0.57295    
Environmentsubsidized          -15.35332    8.57337   69.12686  -1.791  0.07770 .  
Year                             0.50199    0.87435   67.88010   0.574  0.56778    
Sexmale                         17.19067    5.06107   24.45954   3.397  0.00233 ** 
N100                             0.18371    0.04247   52.43447   4.326 6.84e-05 ***
Environmentsubsidized:Sexmale  -12.12232    6.51056   25.59879  -1.862  0.07413 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Envrnm Year   Sexmal N100  
Envrnmntsbs  0.851                            
Year        -1.000 -0.852                     
Sexmale     -0.031  0.298  0.029              
N100         0.069  0.124 -0.070 -0.057       
Envrnmnts:S  0.031 -0.317 -0.030 -0.782  0.126

ANOVA table for 100% MCPs between the two populations

anova(RMmod.year)

Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
Environment      593.69  593.69     1 72.975  6.9058  0.010467 *  
Year              28.34   28.34     1 67.880  0.3296  0.567782    
Sex             1019.81 1019.81     1 24.480 11.8624  0.002073 ** 
N100            1608.73 1608.73     1 52.434 18.7127 6.841e-05 ***
Environment:Sex  298.04  298.04     1 25.599  3.4668  0.074129 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

RM-ANOVA for 95% MCP analyses between the subsidized and non-subsidized

summary(RMmod.year95)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: Home_Range_95mcp ~ Environment + Year + Sex + N95 + Environment *  
    Sex + (1 | Gila)
   Data: year

REML criterion at convergence: 427

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.40965 -0.32439 -0.01753  0.33170  2.10809 

Random effects:
 Groups   Name        Variance Std.Dev.
 Gila     (Intercept) 62.74    7.921   
 Residual             14.71    3.835   
Number of obs: 68, groups:  Gila, 31

Fixed effects:
                                Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                   -847.23191  832.90246   36.03605  -1.017 0.315837    
Environmentsubsidized          -17.57933    5.60239   56.19888  -3.138 0.002710 ** 
Year                             0.43227    0.41610   36.03494   1.039 0.305791    
Sexmale                         18.07841    4.86595   26.11740   3.715 0.000973 ***
N95                              0.03418    0.03342   37.01539   1.023 0.313095    
Environmentsubsidized:Sexmale  -11.52577    5.79029   33.81958  -1.991 0.054663 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Envrnm Year   Sexmal N95   
Envrnmntsbs  0.598                            
Year        -1.000 -0.602                     
Sexmale     -0.025  0.454  0.021              
N95         -0.003  0.254  0.000 -0.027       
Envrnmnts:S  0.004 -0.490 -0.001 -0.841  0.064

ANOVA table for 95% MCPs between the two populations

anova(RMmod.year95)

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
Environment     335.56  335.56     1 61.986 22.8175 1.131e-05 ***
Year             15.87   15.87     1 36.035  1.0792 0.3057907    
Sex             266.80  266.80     1 33.473 18.1420 0.0001568 ***
N95              15.38   15.38     1 37.015  1.0459 0.3130948    
Environment:Sex  58.27   58.27     1 33.820  3.9622 0.0546634 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Raw group means AND adjusted EMMs of Yearly Overall 100%MCP

Directional means of home range (100% MCP) after being adjusted for year, sex and sample size

kable(ref_dfRM_sex, format = "pandoc", caption = 'Table | Directional means of home range (100% MCP) after being adjusted for year, sex and sample size.')

Table | Directional means of home range (100% MCP) after being adjusted for year, sex and sample size.

Environment	Sex	lsmean	SE	df	lower.CL	upper.CL
nonsubsidized	female	23.872887	6.030104	67.38353	11.838009	35.90777
subsidized	female	8.300524	3.305743	48.63495	1.656123	14.94493
nonsubsidized	male	43.544702	6.077155	66.92800	31.414411	55.67499
subsidized	male	13.421268	4.084847	54.24885	5.232502	21.61003

summary(RM.95KDEmod.year)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: Home_Range_95kde ~ Environment + Year + Sex + N + Environment *  
    Sex + (1 | Gila)
   Data: year

REML criterion at convergence: 450.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.48925 -0.42071 -0.07388  0.41738  2.69609 

Random effects:
 Groups   Name        Variance Std.Dev.
 Gila     (Intercept) 175.44   13.25   
 Residual              60.06    7.75   
Number of obs: 61, groups:  Gila, 29

Fixed effects:
                                Estimate Std. Error         df t value Pr(>|t|)   
(Intercept)                   -1.565e+03  1.778e+03  3.477e+01  -0.880  0.38490   
Environmentsubsidized         -2.190e+01  1.132e+01  5.238e+01  -1.934  0.05856 . 
Year                           8.007e-01  8.883e-01  3.479e+01   0.901  0.37356   
Sexmale                        3.244e+01  9.318e+00  2.798e+01   3.481  0.00166 **
N                             -1.727e-02  7.301e-02  3.873e+01  -0.237  0.81428   
Environmentsubsidized:Sexmale -1.805e+01  1.161e+01  2.512e+01  -1.554  0.13264   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Envrnm Year   Sexmal N     
Envrnmntsbs  0.616                            
Year        -1.000 -0.620                     
Sexmale      0.000  0.473 -0.003              
N            0.040  0.357 -0.044  0.074       
Envrnmnts:S  0.017 -0.476 -0.015 -0.796  0.026

anova(RM.95KDEmod.year)

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
Environment     578.08  578.08     1 54.632  9.6250 0.0030363 ** 
Year             48.80   48.80     1 34.786  0.8126 0.3735574    
Sex             956.52  956.52     1 26.005 15.9260 0.0004788 ***
N                 3.36    3.36     1 38.732  0.0559 0.8142796    
Environment:Sex 145.08  145.08     1 25.121  2.4156 0.1326414    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

RM-ANOVA of 95% KDEs for the subsidized population

RM.KDEmod.year<-lmer(Home_Range_95kde~Year+Sex+N+(1|Gila),data = sub)

summary(RM.KDEmod.year)

ANOVA Table for 95%KDE

anova(RM.KDEmod.year)

Type III Analysis of Variance Table with Satterthwaite's method
     Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
Year  33.20   33.20     1 40.016  0.5720 0.453900   
Sex  496.66  496.66     1 17.763  8.5562 0.009132 **
N     21.33   21.33     1 37.608  0.3675 0.548005   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

RM-ANOVA of 50% KDEs for the subsidized

RM.KDE.50.mod.year<-lmer(Home_Range_50kde~Year+Sex+N+(1|Gila),data = sub)

summary(RM.KDE.50.mod.year)

ANOVA Talbe of 50% KDE for the subsidized

anova(RM.KDE.50.mod.year)

TABLE. Raw Group 50% KDE home range means male and female home ranges at Stone Canyon

YR_GRP_Means.50KDE <- summarySE(sub, measurevar="Home_Range_50kde",
                            groupvars=c("Sex"),na.rm = TRUE)

kable(YR_GRP_Means.50KDE, format = "pandoc", caption = 'Table 5 | Raw Group 50% KDE home range means male and female home ranges at Stone Canyon.')

Raw group means AND adjusted EMMs of Yearly Overall 95% KDEs between non/subsidized populations

Collective grid of 100% MCP and 95% KDE of both sites from above

43.4 male 42.9 female Yearly overall means of 95% KDEs grouped by site and sex

YR_Means.95KDEall <- summarySE(year, measurevar="Home_Range_95kde",
                            groupvars=c("Environment","Sex"),na.rm = TRUE)
 
kable(YR_Means.95KDEall, format = "pandoc", caption = 'Table | Raw Group 95% KDE home range means male and female home ranges at non/subsidized.')

Table | Raw Group 95% KDE home range means male and female home ranges at non/subsidized.

Environment	Sex	N	Home_Range_95kde	sd	se	ci
nonsubsidized	female	5	36.80000	9.603905	4.294997	11.924824
nonsubsidized	male	6	69.40000	27.763789	11.334520	29.136310
subsidized	female	37	22.98892	11.046272	1.815996	3.683010
subsidized	male	13	35.00308	12.057546	3.344161	7.286302

Pairwise Comparisons, between sexes by environment, and between environments averaged across sex

emm_s.TK

$emmeans
Sex = female:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  23.76 6.25 68.3    11.29     36.2
 subsidized      8.41 3.47 49.8     1.45     15.4

Sex = male:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  40.95 6.17 69.7    28.65     53.3
 subsidized     13.48 4.27 55.3     4.91     22.0

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Sex = female:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized     15.4 8.72 70.5 1.761   0.0826 

Sex = male:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized     27.5 9.13 71.7 3.009   0.0036 

Graphical Comparisons of Sex Within Each Environment:

plot(emm_s.t2, comparisons = TRUE, xlab = "Least Square Mean (ha)", ylab = "Environment")

Pairwise by sex between enviornements 100% MCP, and 95% KDEs

emm_s.t3

$emmeans
Sex = female:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  23.76 6.25 68.3    11.29     36.2
 subsidized      8.41 3.47 49.8     1.45     15.4

Sex = male:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  40.95 6.17 69.7    28.65     53.3
 subsidized     13.48 4.27 55.3     4.91     22.0

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Sex = female:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized     15.4 8.72 70.5 1.761   0.0826 

Sex = male:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized     27.5 9.13 71.7 3.009   0.0036 

Graphical Comparisons of Sex between the two populations:

plot(emm_s.t3, comparisons = TRUE, xlab = "Least Square Mean (ha)", ylab = "Environment")

Ineractive map of MCPs at Stone Canyon

Create stagnant stamen map of MCPs at Stone Canyon

Error: Don't know how to add north2(SC_stamen_map, x = 0.89, y = 0.85, scale = 0.1, symbol = 16) to a plot

myMap

768x1280 terrain map image from Stamen Maps. 
See ?ggmap to plot it.

Interactive map of KDEs at Stone Canyon

TABLE

Table | Subsidized and non-subsidized directional means of KDE home ranges after being adjusted for year, sex and sample size.

Environment	Sex	lsmean	SE	df	lower.CL	upper.CL
nonsubsidized	female	42.44837	9.346459	53.01148	23.70185	61.19490
subsidized	female	20.93992	4.083131	22.92042	12.49169	29.38814
nonsubsidized	male	80.81873	9.084540	52.50968	62.59348	99.04398
subsidized	male	35.27439	5.501543	29.24014	24.02648	46.52229

SEASONAL ANALYSES

Map of seasonal fluctions of home ranges

## Create MCP polygons by SEASON:
M215_mcp.EM<-mcp_analysis.POLY("./M215/Emergence .csv", percentage= 100)
M215_mcp.DRY<-mcp_analysis.POLY("./M215/Dry .csv", percentage= 100)
M215_mcp.MON<-mcp_analysis.POLY("./M215/Monsoon .csv", percentage= 100)

M112_mcp.DRY<-mcp_analysis.POLY("./M112/Dry .csv", percentage= 100)
M112_mcp.MON<-mcp_analysis.POLY("./M112/Monsoon .csv", percentage= 100)
M112_mcp.PM<-mcp_analysis.POLY("./M112/Post_Monsoon .csv", percentage= 100)

M119_mcp.DRY<-mcp_analysis.POLY("./M119/Dry .csv", percentage= 100)
M119_mcp.MON<-mcp_analysis.POLY("./M119/Monsoon .csv", percentage= 100)
M119_mcp.PM<-mcp_analysis.POLY("./M119/Post_Monsoon .csv", percentage= 100)

F114_mcp.EM<-mcp_analysis.POLY("./F114/Emergence .csv", percentage= 100)
F114_mcp.DRY<-mcp_analysis.POLY("./F114/Dry .csv", percentage= 100)
F114_mcp.MON<-mcp_analysis.POLY("./F114/Monsoon .csv", percentage= 100)
F114_mcp.PM<-mcp_analysis.POLY("./F114/Post_Monsoon .csv", percentage= 100)

F137_mcp.EM<-mcp_analysis.POLY("./F137/Emergence .csv", percentage= 100)
F137_mcp.DRY<-mcp_analysis.POLY("./F137/Dry .csv", percentage= 100)
F137_mcp.MON<-mcp_analysis.POLY("./F137/Monsoon .csv", percentage= 100)
F137_mcp.PM<-mcp_analysis.POLY("./F137/Post_Monsoon .csv", percentage= 100)

F147_mcp.EM<-mcp_analysis.POLY("./F147/Emergence .csv", percentage= 100)
F147_mcp.DRY<-mcp_analysis.POLY("./F147/Dry .csv", percentage= 100)
F147_mcp.MON<-mcp_analysis.POLY("./F147/Monsoon .csv", percentage= 100)
F147_mcp.PM<-mcp_analysis.POLY("./F147/Post_Monsoon .csv", percentage= 100)

F252_mcp.EM<-mcp_analysis.POLY("./F252/Emergence .csv", percentage= 100)
F252_mcp.DRY<-mcp_analysis.POLY("./F252/Dry .csv", percentage= 100)
F252_mcp.MON<-mcp_analysis.POLY("./F252/Monsoon .csv", percentage= 100)
F252_mcp.PM<-mcp_analysis.POLY("./F252/Post_Monsoon .csv", percentage= 100)

F36_mcp.EM<-mcp_analysis.POLY("./F36/Emergence .csv", percentage= 100)
F36_mcp.DRY<-mcp_analysis.POLY("./F36/Dry .csv", percentage= 100)
F36_mcp.MON<-mcp_analysis.POLY("./F36/Monsoon .csv", percentage= 100)
F36_mcp.PM<-mcp_analysis.POLY("./F36/Post_Monsoon .csv", percentage= 100)

F66_mcp.EM<-mcp_analysis.POLY("./F66/Emergence .csv", percentage= 100)
F66_mcp.DRY<-mcp_analysis.POLY("./F66/Dry .csv", percentage= 100)
F66_mcp.MON<-mcp_analysis.POLY("./F66/Monsoon .csv", percentage= 100)
F66_mcp.PM<-mcp_analysis.POLY("./F66/Post_Monsoon .csv", percentage= 100)

## Fortify mcp polygons for ggplot2 *SEASON*:
M215_mcp.EMT <- fortify(M215_mcp.EM, region = "id")
M215_mcp.DRYT <- fortify(M215_mcp.DRY, region = "id")
M215_mcp.MONT <- fortify(M215_mcp.MON, region = "id")

M112_mcp.DRYT <- fortify(M112_mcp.DRY, region = "id")
M112_mcp.MONT <- fortify(M112_mcp.MON, region = "id")
M112_mcp.PMT <- fortify(M112_mcp.PM, region = "id")

M119_mcp.DRYT <- fortify(M119_mcp.DRY, region = "id")
M119_mcp.MONT <- fortify(M119_mcp.MON, region = "id")
M119_mcp.PMT <- fortify(M119_mcp.PM, region = "id")

F114_mcp.EMT <- fortify(F114_mcp.EM, region = "id")
F114_mcp.DRYT <- fortify(F114_mcp.DRY, region = "id")
F114_mcp.MONT <- fortify(F114_mcp.MON, region = "id")
F114_mcp.PMT <- fortify(F114_mcp.PM, region = "id")

F137_mcp.EMT <- fortify(F137_mcp.EM, region = "id")
F137_mcp.DRYT <- fortify(F137_mcp.DRY, region = "id")
F137_mcp.MONT <- fortify(F137_mcp.MON, region = "id")
F137_mcp.PMT <- fortify(F137_mcp.PM, region = "id")

F147_mcp.EMT <- fortify(F147_mcp.EM, region = "id")
F147_mcp.DRYT <- fortify(F147_mcp.DRY, region = "id")
F147_mcp.MONT <- fortify(F147_mcp.MON, region = "id")
F147_mcp.PMT <- fortify(F147_mcp.PM, region = "id")

F252_mcp.EMT <- fortify(F252_mcp.EM, region = "id")
F252_mcp.DRYT <- fortify(F252_mcp.DRY, region = "id")
F252_mcp.MONT <- fortify(F252_mcp.MON, region = "id")
F252_mcp.PMT <- fortify(F252_mcp.PM, region = "id")

F36_mcp.EMT <- fortify(F36_mcp.EM, region = "id")
F36_mcp.DRYT <- fortify(F36_mcp.DRY, region = "id")
F36_mcp.MONT <- fortify(F36_mcp.MON, region = "id")
F36_mcp.PMT <- fortify(F36_mcp.PM, region = "id")

F66_mcp.EMT <- fortify(F66_mcp.EM, region = "id")
F66_mcp.DRYT <- fortify(F66_mcp.DRY, region = "id")
F66_mcp.MONT <- fortify(F66_mcp.MON, region = "id")
F66_mcp.PMT <- fortify(F66_mcp.PM, region = "id")

mcp.shift.TEST5 <- ggplot() +
  geom_polygon(data=F114_mcp.EMT, aes(x=F114_mcp.EMT$long, y=F114_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F114_mcp.DRYT, aes(x=F114_mcp.DRYT$long, y=F114_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F114_mcp.MONT, aes(x=F114_mcp.MONT$long, y=F114_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F114_mcp.PMT, aes(x=F114_mcp.PMT$long, y=F114_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F137_mcp.EMT, aes(x=F137_mcp.EMT$long, y=F137_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F137_mcp.DRYT, aes(x=F137_mcp.DRYT$long, y=F137_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F137_mcp.MONT, aes(x=F137_mcp.MONT$long, y=F137_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F137_mcp.PMT, aes(x=F137_mcp.PMT$long, y=F137_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F147_mcp.EMT, aes(x=F147_mcp.EMT$long, y=F147_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F147_mcp.DRYT, aes(x=F147_mcp.DRYT$long, y=F147_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F147_mcp.MONT, aes(x=F147_mcp.MONT$long, y=F147_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F147_mcp.PMT, aes(x=F147_mcp.PMT$long, y=F147_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  # geom_polygon(data=F252_mcp.EMT, aes(x=F252_mcp.EMT$long, y=F252_mcp.EMT$lat),
  #              alpha=0.1,colour="black",linetype=2) +
  # geom_polygon(data=F252_mcp.DRYT, aes(x=F252_mcp.DRYT$long, y=F252_mcp.DRYT$lat),
  #              alpha=0.1,colour="black",linetype=3) +
  # geom_polygon(data=F252_mcp.MONT, aes(x=F252_mcp.MONT$long, y=F252_mcp.MONT$lat),
  #              alpha=0.1,colour="black",linetype=4) +
  # geom_polygon(data=F252_mcp.PMT, aes(x=F252_mcp.PMT$long, y=F252_mcp.PMT$lat),
  #              alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F36_mcp.EMT, aes(x=F36_mcp.EMT$long, y=F36_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F36_mcp.DRYT, aes(x=F36_mcp.DRYT$long, y=F36_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F36_mcp.MONT, aes(x=F36_mcp.MONT$long, y=F36_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F36_mcp.PMT, aes(x=F36_mcp.PMT$long, y=F36_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  geom_polygon(data=F66_mcp.EMT, aes(x=F66_mcp.EMT$long, y=F66_mcp.EMT$lat),
               alpha=0.1,colour="blue",linetype=2) +
  geom_polygon(data=F66_mcp.DRYT, aes(x=F66_mcp.DRYT$long, y=F66_mcp.DRYT$lat),
               alpha=0.1,colour="red",linetype=3) +
  geom_polygon(data=F66_mcp.MONT, aes(x=F66_mcp.MONT$long, y=F66_mcp.MONT$lat),
               alpha=0.1,colour="green",linetype=4) +
  geom_polygon(data=F66_mcp.PMT, aes(x=F66_mcp.PMT$long, y=F66_mcp.PMT$lat),
               alpha=0.1,colour="black",linetype=5) +
  theme_bw() +
  labs(x="Easting (m)", y="Northing (m)") +
  labs(caption = "Figure 6 |  SC seasonal home range shifts of five lizards. All seasonal polygons stay relatively stable with \n considerable overlap and without any major shifts.")+
  theme(plot.caption = element_text(hjust = 0,lineheight = 0.9))+
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))

mcp.shift.TEST5

TABLE group means of seasonal home ranges between the two populations averaged across sex

seasonal<-read.csv("SC_Seasonal_Data.csv")

library(Rmisc)

SEAS_GRP_Means <- summarySE(seasonal, measurevar="Home_Range_100mcp",
                            groupvars=c("Environment","Season"), na.rm = TRUE)

# SEAS_GRP_Means
kable(SEAS_GRP_Means, format = "pandoc", caption = 'Table 6 | Group means of seasonal home ranges between Stone Canyon (subsidized) and Owl Head Buttes (non-subsidized). These means are averaged across sex.')

RM-ANOVA for seasonal home ranges between environments

library(lme4)
library(readr)
library(lmerTest)
# seasonal<-read.csv("SC_Seasonal_Data.csv")

RM.mod.Season <- lmer(Home_Range_100mcp~Environment+Season+Sex+N+Environment*Season+(1|Gila), 
                      data=seasonal)
summary(RM.mod.Season)

ANOVA table of seasonal HRs between envs.

anova(RM.mod.Season)

TABLE of seasonal home ranges by sex between the two populations

SEAS_GRP_TEST <- summarySE(seasonal, measurevar="Home_Range_100mcp",
                           groupvars=c("Environment","Season","Sex"), na.rm = TRUE)
# SEAS_GRP_Means
kable(SEAS_GRP_TEST, format = "pandoc", caption = 'Table 7 | Seasonal home range means between Stone Canyon (subsidized) and Owl Head Buttes (non-subsidized) popuations for males and females. These are raw means before being adjusted for environment, season, sex, and sample size.')

Table 7 | Seasonal home range means between Stone Canyon (subsidized) and Owl Head Buttes (non-subsidized) popuations for males and females. These are raw means before being adjusted for environment, season, sex, and sample size.

Environment	Season	Sex	N	Home_Range_100mcp	sd	se	ci
nonsubsidized	Dry	female	5	15.6600000	8.6291946	3.8590932	10.7145603
nonsubsidized	Dry	male	7	29.4714286	12.6476235	4.7803524	11.6971008
nonsubsidized	Emergence	female	5	4.4600000	3.4333657	1.5354478	4.2630866
nonsubsidized	Emergence	male	5	1.1600000	1.8242807	0.8158431	2.2651436
nonsubsidized	Monsoon	female	6	22.9833333	9.8151753	4.0070285	10.3003948
nonsubsidized	Monsoon	male	7	24.2285714	9.8668999	3.7293376	9.1253605
nonsubsidized	Post_Monsoon	female	4	1.4000000	1.4491377	0.7245688	2.3059014
nonsubsidized	Post_Monsoon	male	7	0.2857143	0.3670993	0.1387505	0.3395102
subsidized	Dry	female	11	10.1754545	8.0883118	2.4387178	5.4338018
subsidized	Dry	male	6	18.2816667	13.2661214	5.4158714	13.9219406
subsidized	Emergence	female	6	2.1133333	1.8474920	0.7542354	1.9388239
subsidized	Emergence	male	3	2.0666667	1.5326556	0.8848792	3.8073277
subsidized	Monsoon	female	11	10.6918182	8.4988679	2.5625051	5.7096172
subsidized	Monsoon	male	7	10.3528571	6.3010018	2.3815548	5.8274547
subsidized	Post_Monsoon	female	11	3.6309091	5.5527983	1.6742317	3.7304207
subsidized	Post_Monsoon	male	3	0.6333333	0.8007705	0.4623250	1.9892241

figures for raw seasonal home ranges between the two populations

Figures Adjusted EMMs of seasonal home range between the two populations

Collective grid of raw and adjusted seasonal home ranges

ggarrange(raw.seasonal, adj.seasonal, labels = c("A", "B"),
          nrow = 2)

Post hoc analyses of seasonal home ranges

Pairwise of each season between populations, overaged over levels of sex

emm_s.t <- emmeans(RM.mod.Season, pairwise ~ Environment | Season)
emm_s.t

$emmeans
Season = Dry:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  18.66 2.13 85.1   14.428    22.89
 subsidized     13.21 1.68 84.3    9.866    16.56

Season = Emergence:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized   3.36 2.11 85.3   -0.844     7.56
 subsidized      3.99 2.44 83.5   -0.864     8.84

Season = Monsoon:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized  21.77 1.92 83.7   17.943    25.59
 subsidized      8.60 1.69 82.8    5.239    11.96

Season = Post_Monsoon:
 Environment   emmean   SE   df lower.CL upper.CL
 nonsubsidized   3.21 2.26 83.7   -1.290     7.72
 subsidized      3.99 2.08 84.1   -0.152     8.13

Results are averaged over the levels of: Sex 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Season = Dry:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized    5.448 2.67 85.0  2.037  0.0448 

Season = Emergence:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized   -0.630 3.16 85.1 -0.199  0.8426 

Season = Monsoon:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized   13.167 2.52 83.5  5.219  <.0001 

Season = Post_Monsoon:
 contrast                   estimate   SE   df t.ratio p.value
 nonsubsidized - subsidized   -0.776 2.99 86.4 -0.259  0.7960 

Results are averaged over the levels of: Sex 

Graphical comparisons

plot(emm_s.t, comparisons = TRUE)

Pairwise between seasons within each popultion

emm_s.t4

$emmeans
Environment = nonsubsidized:
 Season       emmean   SE   df lower.CL upper.CL
 Dry           18.86 2.25 88.4   14.383    23.34
 Emergence      3.33 2.24 88.7   -1.118     7.77
 Monsoon       21.85 2.03 87.5   17.811    25.89
 Post_Monsoon   2.36 2.36 87.0   -2.322     7.04

Environment = subsidized:
 Season       emmean   SE   df lower.CL upper.CL
 Dry           12.23 1.75 87.4    8.745    15.72
 Emergence      4.32 2.55 84.7   -0.741     9.39
 Monsoon        9.04 1.78 86.0    5.515    12.57
 Post_Monsoon   5.09 2.07 85.8    0.981     9.21

Results are averaged over the levels of: Sex 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Environment = nonsubsidized:
 contrast                 estimate   SE   df t.ratio p.value
 Dry - Emergence            15.532 3.07 69.4  5.054  <.0001 
 Dry - Monsoon              -2.992 2.89 67.3 -1.036  0.7292 
 Dry - Post_Monsoon         16.500 3.24 78.9  5.098  <.0001 
 Emergence - Monsoon       -18.524 2.91 68.0 -6.361  <.0001 
 Emergence - Post_Monsoon    0.968 3.08 73.0  0.314  0.9891 
 Monsoon - Post_Monsoon     19.492 3.03 74.0  6.426  <.0001 

Environment = subsidized:
 contrast                 estimate   SE   df t.ratio p.value
 Dry - Emergence             7.907 3.11 88.6  2.543  0.0602 
 Dry - Monsoon               3.187 2.28 66.0  1.395  0.5070 
 Dry - Post_Monsoon          7.137 2.68 80.2  2.666  0.0450 
 Emergence - Monsoon        -4.720 3.20 89.6 -1.475  0.4569 
 Emergence - Post_Monsoon   -0.769 2.94 77.2 -0.262  0.9937 
 Monsoon - Post_Monsoon      3.951 2.78 84.9  1.421  0.4899 

Results are averaged over the levels of: Sex 
P value adjustment: tukey method for comparing a family of 4 estimates 

Graphical Comps

plot(emm_s.t4, comparisons = TRUE)

Pairwise between sexes of each season of the subsidized population

sub <- subset(seasonal, Environment == "subsidized")

RM.mod.Sub <- lmer(Home_Range_100mcp~Season+Sex+N+Season*Sex+(1|Gila), data=sub)

emm_s.t5 <- emmeans(RM.mod.Sub, pairwise ~ Sex | Season)
emm_s.t5 

Graphical Comps

plot(emm_s.t5, comparisons = TRUE)

Pairwise between sexes of each season of the non-subsidized population

nonsub <- subset(seasonal, Environment == "nonsubsidized")
View(nonsub)
RM.mod.NSub <- lmer(Home_Range_100mcp~Season+Sex+N+Season*Sex+(1|Gila), data=nonsub)

emm_s.t6 <- emmeans(RM.mod.NSub, pairwise ~ Sex | Season)
emm_s.t6 

Graphical Comps

plot(emm_s.t6, comparisons = TRUE)